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^^ I We calculate numerically the spectrum of ultra-high energy cosmic rays on Earth 

^ ■ assuming that their sources are distributed in space like the observed galaxies. We 

^«0 , use the CfA2 and the PSCz galaxy redshift surveys to model the local galaxy distri- 

^1^ \ bution, properly taking into account the galaxy selection functions for each survey. 

Q^ ■ When the survey selection effects are included, we find that the local overdensity is 

^^ , only a factor of two, an order of magnitude less than used in some earlier studies. An 

Q \ overdensity of two is not enough to bridge the gap between the predicted number 

of events above 10^° eV and the measured flux at these highest energies. This con- 

O . elusion is particularly strong for soft injection spectra (ex E^^) where the observed 

JL \ number of events is 6 u higher than the expected one. However, if the injection 

$-H ' spectrum is hard (oc E~'^), the small local overdensity helps bring the present data 

c/3 . within 2a of the low number of events predicted above 10^'^ eV. In this case, the 

. . i Greisen-Zatzepin-Kuzmin cutoff is not a cutoff but rather a feature in the cosmic 

.5^ ' ray spectrum. 
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1 Introduction 



The unexpected detection of cosmic rays with energies above 10^° eV has 
triggered considerable interest in the possible origin and nature of these par- 
ticles [1] . These highest energy events are surprising for the following reasons. 
If these particles are protons, they likely originate in extragalactic sources, 
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since at these high energies the Galactic magnetic field cannot confine pro- 
tons in the Galaxy. However, extragalactic protons with energies above a few 
times 10^^ eV can produce pions through interactions with the cosmic mi- 
crowave background (CMB) and consequently lose significant amounts of en- 
ergy as they traverse intergalactic distances [2]. Therefore, in addition to the 
extraordinary energy requirements for astrophysical sources to accelerate pro- 
tons to ^10^° eV, the photopion threshold reaction suppresses the observable 
fiux above ~ 10^'^ eV. These conditions were expected to cause a natural 
high-energy limit to the cosmic ray spectrum known as the Greisen-Zatsepin- 
Kuzmin (GZK) cutoff [2]. 

As clearly shown by the most recent compilation of the AGASA collaboration 
data [3] , the spectrum of cosmic rays does not end at the expected GZK cutoff. 
The significant fiux observed above 10^° eV together with a nearly isotropic 
distribution of event arrival directions has challenged astrophysically based 
scenarios (see e.g., [4] and references therein) and has inspired a number of 
alternatives (see, e.g., [5] and references therein). However, the GZK cutoff 
is not an absolute end to the cosmic ray spectrum but it generates a feature 
around 5 x 10^^ eV. The spectrum recovers at energies above this feature [6] 
and the local distribution of sources can significantly affect the agreement 
between predicted and observed spectra. In particular, a local overdensity of 
sources can decrease the gap between observed and detected events above the 
GZK cutoff [7-9] . This effect can easily be understood: photopion energy losses 
limit the maximum distance at which sources can contribute at the highest 
energies to a few tens of Mpc, while cosmic rays below the pion production 
threshold come from much larger volumes. A local overdensity will increase 
the observed flux at the highest energies relative to the lower energy flux. Here 
we show that the observed local overdensity is not high enough to explain the 
data unless, perhaps, the sources have a hard injection spectrum. 

In the following, we simulate the propagation of ultra-high energy cosmic rays 
(UHECRs) from extragalactic sources to Earth. Our numerical propagation 
code includes pair production, photopion production, and adiabatic losses due 
to the universe's expansion. For a uniform distribution of sources, our numer- 
ical results agree well with previous results and analytical calculations. We 
study the effect of a local overdensity by assuming that the number density of 
sources is proportional to the number density of galaxies. To keep this study 
independent of speciflc hypothetical sources, we use different injection spec- 
tra and consider a possible source evolution with redshift. However, since we 
assume that UHECRs are protons, our results are mostly relevant for extra- 
galactic astrophysical acceleration scenarios. 

Rather than adopting analytic models of the galaxy distribution, as do [9,10], 
we extract the distribution of galaxies from observations of large scale struc- 
ture using the CfA2 and the PSCz galaxy catalogs. We assume that the den- 



sity field derived from tliese studies lias the same shape as the density field 
of UHECR sources, although a priori UHECR sources may cluster differently 
from luminous matter. We show that the local density is only about a factor 
of two above the mean, in contrast with much higher estimates previously 
published [8]. This large discrepancy was caused by neglecting the necessary 
galaxy selection functions which account for the fact that nearby galaxies are 
far easier to detect than far away galaxies. Once we include the selection func- 
tions, we find that the real overdensity is not high enough to bridge the gap 
between predicted and observed spectra for soft sources {J{E) oc E~'^, with 
7 = 3). However, sources with a hard injection spectrum (7^2.1) can fit the 
present data within 2a, at energies above a few 10^^ eV. Sources of UHECRs 
distributed as ordinary galaxies are marginally consistent with present data 
and, for hard injection spectra, the GZK cutoff is not really a cutoff but a 
feature in the high-energy cosmic ray spectrum. 

The plan of this paper is as follows. In §2, we discuss the proper way to model a 
distribution of sources associated with galaxies. We derive the galaxy density 
field in our neighborhood of the universe using the CfA2 and PSCz galaxy 
redshift surveys with their respective selection function. In §3, we describe our 
UHECR propagation code. In §4, we display the results for different injection 
spectra and a realistic spatial distribution of sources. We also contrast our 
results with previous work and analytical estimates. We conclude in §5. 



2 The Galaxy Density Field 



Although galaxies are not homogeneously distributed in the local universe, 
the GZK cutoff is usually calculated assuming a homogeneous distribution 
of sources throughout space. If the sources are distributed like the luminous 
matter around us, the effects of an inhomogeneous galaxy distribution needs 
to be taken into account when predicting the spectrum of UHECRs. In order 
to include the effects of the spatial inhomogeneity in the UHECR spectrum, 
we must consider the estimated galaxy density field in our neighborhood of 
the universe. The density field is usually measured by selecting galaxies from 
an imaging survey of the sky and taking their redshifts. Almost invariably, the 
galaxies are selected to be brighter than some limiting fiux in some band, /um, 
expressed as an "apparent magnitude," mum = — 2.51og]^Q(/iim//o), where /o is 
an arbitrary zero-point. For all (or for some random subsample) of the galaxies 
brighter than this, one takes their spectra and determines their redshifts z. 
According to the Hubble law, the redshifts are related to their distances d = 
Hqcz, where c is the speed of light and Hq = lOOh Mpc/ km/s, with h ~ 0.5- 
0.8 (e.g., [11]). 

However, a fiux-limited survey is not a volume-limited survey. In a fiux-limited 



survey, the raw distribution of redshifts cannot be used without regard to 
the way galaxies were selected. Here we describe the proper way to derive 
density fields from galaxy redshift surveys. We limit ourselves to measuring 
the density in redshift space and do not include the small effects of deviations 
from the Hubble law due to galaxy peculiar velocities. We base our approach 
on methods dating back to [12] (see, e.g., [13] for a recent review). 

We will use two surveys for our analysis. To compare directly with [8], we 
consider the Center for Astrophysics Redshift Survey (CfA2; [14]). Although 
this survey comprises about 10,000 galaxy redshifts (selected to be brighter 
than m = 15.5, a 5-band magnitude), it covers only about 17% of the whole 
sky. In order to better evaluate the effects of the density field of galaxies on 
the cosmic ray spectrum, we should use surveys which probe the density field 
over nearly the whole sky. The best sample of galaxies to use for this purpose 
is the IRAS PSC2; Survey [15], which consists of about 15,000 galaxies with 
infrared fiuxes > 0.6 Jy and covers about 84% of the sky. 

A consequence of the fiux limits in any survey is that at different redshifts, a 
different set of galaxy luminosities L is observed, determined by the faintest 
luminosity observable at that redshift L^-ai^^z). For an Euclidean metric, this 
luminosity is related to the fiux limit by: 



A7r{Hocz)^fi 



lim- 



At cosmological distances (generally only appropriate when z > 0.1) more 
complicated relations apply, which in general depend on the cosmological 
model [16] (see [17] for a useful compilation of results). If the distribution 
of galaxy luminosities is described by the galaxy luminosity function ^{L), 
then the fraction of all galaxies that are observable at any redshift is given by: 



lo^dL^L) ■ 



(2) 



The quantity (f){z) is usually referred to in the literature as the "selection func- 
tion" [18]. The most common methods used to determine the galaxy luminosity 
function from the survey itself are those of [19] and [20]. These methods as- 
sume that the luminosity function is universal (i.e., independent of redshift) 
and use nearby galaxies to determine the shape of the faint end and far away 
galaxies for the shape of the bright end. 

Consider, for example. Figure 1, which shows in the top panel the distribution 
of galaxies and redshifts in CfA2. Here we express galaxy luminosity in terms 
of the "absolute magnitude," M = —2.5 logj^Q L+ const; thus, in the figure, the 
faintest galaxies are at the top. The thick solid line shows the fiux limit of the 



survey, translated into an absolute magnitude limit at each redshift. Because 
of this limit a number of galaxies which are observable at low redshifts are too 
faint to be observed at higher redshift. The fraction of galaxies (f){z) between 
absolute magnitudes —22 < M < —10 which are unobservable at each redshift 
is shown in the bottom panel of Figure 1, based on a fit to the luminosity 
function in the survey using the method of [19]. Because the function falls 
rapidly from unity, it is clear that even at low redshifts the effects of the flux 
limit are important. 
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Fig. 1. The top panel shows the absolute magnitudes (related to luminosity by 
M = — 2.51ogiQ-L+ const) and redshifts of CfA2 galaxies. Shown as the thick solid 
line is the flux limit, converted to the appropriate absolute magnitude at each red- 
shift. The bottom panel shows the fraction of galaxies in the range —22 < M < —10 
that we estimate to be brighter than the flux limit. This function falls rapidly with 
redshift. When interpreting the top plot, remember that the volume probed at low 
redshift is far smaller than that probed at high redshift. 

In principle, we can construct a volume-limited sample by choosing only galax- 
ies brighter than some magnitude M up to the redshift at which the thick solid 
line representing the flux limit crosses M . This usually decreases significantly 
the number of galaxies and the redshift depth of the survey. A more effec- 
tive approach is to use 0(z) to calculate the distribution of observed galaxies 
with redshift which one would expect if the distribution were homogeneous. 
Then, the density field can be inferred by comparing the actual counts to these 



expected counts. The top panel of Figure 2 compares these expected counts 
(dotted hne) in redshift shells of thickness 0.001 to the observed counts (solid 
line) in CfA2. It appears that locally we are in an overdensity of galaxies of 
about a factor of two; note that at large distances, where each shell corre- 
sponds to a considerable amount of volume, and thus averages over large-scale 
structures, the number of galaxies is very nearly the expected number. Instead, 
if the flux limit is neglected (dashed line), the "expected" number of galaxies 
in each shell scales approximately as the square of the redshift of that shell. If 
we normalize these "expected" counts in approximately the same way as [8], 
we recover the incorrect result that we live in an overdensity of approximately 
a factor of 30. It is clear from the comparison of the dotted to the dashed 
curve that neglecting the flux limits is a poor approximation. Put simply, we 
do not live in a large overdensity, but we can certainly detect galaxies more 
easily if they are close by rather than far away. 






100 ^ 



10 ^ 






100 ^ 



10 ^ 



CfA2 Survey 

Huchra, Geller, &: Corwin (1995). 

(17% of the sky) 




f-f 



- PSCz Survey 

Saunders et aL (2000) 
(84% of the sky~ 



Actual Counts 
Expected Counts 
(Based on Flux Limits) 
lixpected Counts 
(Neglecting Flux Limits) 




Fig. 2. CoTnparison of observed counts (solid line) to those predicted based on the 
flux limits (dotted line) and those predicted neglecting the flux limits (dashed line). 
The CfA2 survey is shown at top, the PSCz at bottom. Both show a local overdensity 
of only about a factor of two when the flux limits are properly accounted for. 



As mentioned above, the CfA2 survey covers a relatively small fraction of the 
sky. Thus, the PSCz redshift survey provides a more useful sample to use in 
the context of this paper. Using the selection function provided by [15], we 



again show the expected versus the observed counts for the PSC2; survey in 
the bottom panel of Figure 2. This survey also shows we are living in a slight 
overdensity, and furthermore reveals the general homogeneity of the nearby 
universe. (The actual counts and their dependence on redshift are different 
from CfA2, because the galaxies are selected in different ways). 



3 The Propagation Code for UHECRs 



Armed with a more realistic model of the local universe, we can calculate 
the spectrum of UHECRs that would be observed on Earth for extragalactic 
sources with different injection spectra. Our numerical propagation code in- 
cludes pair production and photopion production as energy losses, as well as 
adiabatic losses due to the expansion of the universe. We compare our numeri- 
cal results with analytical calculations which we carry out as in [6]. In order to 
isolate the effect of density inhomogeneities we neglect the effect of magnetic 
fields in this paper. 

We compare the results of the simulations with the observed spectrum by 
requiring that the total number of simulated events with energy above a nor- 
malizing energy, Enorm, be the same as what is observed. In general we choose 
Enorm — 10^^ cV siucc at lower energies the flux of cosmic rays is likely to 
have a galactic origin and at higher energies the observed number of events 
is very small. For a given source spectral index, we generate many realiza- 
tions of the spectrum on Earth. Each realization has the same total number 
of observed particles above Enorm calculated as follows. The flux of UHE- 
CRs contributed from a shell of thickness A2; at redshift z is proportional 
to p{z) oc (1 + z)™^^/^/(2;), where the redshift dependence of the density of 
sources and the flux suppression due to distance are included explicitly [21]. In 
this formula, we allow for source luminosity evolution through the parameter 
m. The function f{z) describes possible deviations from a homogeneous spatial 
distribution of sources. For a homogeneous distribution f{z) = 1, otherwise 
f{z) represents the local overdensity of sources at redshift z as, for example, 
those derived from the catalogs introduced in §2. We assign the redshift of one 
event by extracting a random number according to the distribution p{z) given 
above. The energy of the particle at the source is extracted from a distribution 
representing the spectrum of the source assumed to have the form of a power 
law E^'^ . The particle is then propagated from the source to the detector. 

The photopion energy losses are simulated following [22]. In each spatial step 
of the simulation of size As (~ 200 kpc), the average number of photons that 
can induce the production of pions in the scattering with a proton of energy 



E at time t is 



<'^''(^-^^")-K^- P' 



Here, 1{E) = c[{l/E){dE/dt)] is taken from [6] and Kp{E) is an effective 
inelasticity at energy E, which can be approximated by 



with Eth = 3 X fO^^ GeV [22]. Once the average number of pion producing 
photons has been determined over the path As, the actual number of photons 
with which the proton interacts is extracted from a Poisson distribution with 
mean {Nph{E, As)). For sufficiently small As, the number of photons in each 
step is usually either zero or one. The energy of each photon is extracted from a 
Planck distribution at temperature Tqmb = 2.728 K (with a minimum energy 
corresponding to the threshold for pion production). Since the CMB photons 
are isotropically distributed in space, we extract the interaction angle from a 
flat distribution in cos(^). The final energy of the proton after each scattering 
with a CMB photon is calculated from the kinematics of the scattering. Since 
the inelasticity for pair production is very low, we treat it as a continuous 
energy loss process. 

We compare our numerical results with analytic results for the modification 
factor from single sources and from a distribution of sources as in [6]. The 
agreement is excellent and the effect of the fiuctuations at energies larger than 
~ (3 — 4) X 10^^ eV is evident in Figures 5-8 as we discuss below. The average 
of the simulated fiux is slightly larger than the analytic one, as expected for 
the stochastic process of photopion production (on small distances there is an 
appreciable chance that some protons do not interact at all). 



4 Results 



As we stressed above, a local overdensity increases the observed fiux above 
GZK energies with respect to the fiux at lower energies. To show this effect 
clearly, we estimated the change in the spectrum due to a simple top-hat 
model before considering more realistic models of the galaxy density field. In 
Figure 3, we display the results of our analytical calculation for three choices 
of the overdensity, p/p = 1, 10, 30 (solid, dashed, and dash-dotted lines 
respectively), all in a volume of radius ~ 20 Mpc around the Earth and with 
source spectral index 7 = 3. From this figure, it is clear that the overdensity 



increases the flux at the highest energies versus the flux at and below the GZK 
feature, as mentioned above. 
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Fig. 3. The effect of a local overdensity within 20 Mpc on the fluxes of UHE- 
CRs, according to our analytical calculation for three choices of the overdensity, 
pIp = 1, 10, 30 (solid, dashed and dash-dotted lines respectively), all in a volume 
of radius ~ 20 Mpc around the Earth and with source spectral index 7 = 3. 

The galaxy surveys discussed in §2 provide a more realistic model of the local 
density field. Among the different surveys, the PSCz catalog covers a lot more 
solid angle and reaches further (up to a redshift z^ax = 0.1) when compared 
to the CfA2 survey. Thus, we use the PSCz to study the effect of different 
source spectra below. But first we compare the results of the two catalogs for 
a fixed spectral index (7 = 3) with the homogeneous distribution in Figure 
4. This choice of 7 allows a direct comparison of our results with those of 
[8]. In this figure, we normalize the spectrum by requiring that the number of 
events with energy above 10^^ eV equal the AGASA number of 728 [23,24]. 
The error bars in the simulation are obtained by generating 100 realizations 
and calculating the mean and variance of the set. Figure 4 shows that the two 
catalogs give very similar results. 

In Figure 5, we compare the case of a homogeneous distribution and the PSCz 
catalog versus the AGASA data also for 7 = 3. The total number of events 
with E > 10^° eV is 1.2 ± 1.0 for the homogeneous case and 1.5 ± 1.0 for the 
PSCz catalog of galaxies. Both numbers are consistent with 1, while AGASA 
has detected 8 events with E > 10^'^ eV. The data is more than 6 a away from 
the observations. None of our realizations have the observed number of events 
for this soft spectrum. In the figure, the solid and dashed lines represent the 
result of the analytical calculations for the same value of the parameters and 
for the homogeneous and PSCz cases respectively. The dash-dotted and dash- 




Fig. 4. Simulated fluxes of UHECRs for 7 = 3 using the CfA2 catalog (hatched 
upwards to the right, i.e., like "/") and the PSCz catalog (hatched downwards to the 
right, i.e., like '\"). 

dot-dot-dotted lines trace the mean simulated fluxes for the homogeneous case 
(hatched downwards to the right, i.e., like "\") and the PSCz case (hatched 
upwards to the right, i.e., like "/"). As can be seen from Figure 5, the difference 
between the analytic calculations and the mean of the numerical calculation is 
quite small. In addition, the homogeneous case has only slightly lower fluxes at 
higher energies than the realistic models. Here our results differ significantly 
from those of [8] , due to the fact that the real local overdensity is much smaller 
than found in [8]. 

In Figure 6, we show the generated fluxes assuming that 9075 events have 
been observed above 10^^ eV, a number which the Auger project is expected 
to reach in the first 3 years of operation. Again we compare the results of a 
homogeneous distribution of sources (\ hatches) with those associated with 
the PSCz galaxy distribution (/ hatches), using 7 = 3 for both. As expected 
from a larger sample, the size of the error bars decreases significantly. The 
number of events at E > 10^'^ eV in the two cases is still much smaller than if 
extrapolated from the AGASA observations (~ 100). This figure shows that 
the Auger project will be able to constrain the nature of the GZK feature 
much more accurately than possible at present. 

It is clear that with a source injection spectrum as steep as 7 = 3, the number 
of events with E > 10^*^ eV is significantly smaller than what has been observed 
so far. Since sources with harder spectra have been proposed, we consider 
different choices for the injection spectrum. We simulated the case of 7 = 2.7, 
but the results do not diffr significantly from the 7 = 3 case. Next we discuss 
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Fig. 5. Simulated fluxes for the AG AS A statistics of 728 events above 10"^^ eV, 
and 7 = 3, using a homogeneous source distribution (\ hatches) and the PSCz 
distribution (/ hatches). The solid and dashed lines are the results of the analytical 
calculations for the same two cases. The dash-dotted and dash- dot- dot- dotted lines 
trace the mean simulated fluxes for the homogeneous and the PSCz cases. 

the case of 7 = 2.1, where we find the gap between observed and predicted 
fiux to reach the 2 a level. 

In Figure 7, we show the results for 7 = 2.1 for three cases: a homogeneous 
distribution of the sources with Zmax = 0.1 (/ hatches), the PSCz distribution 
with Zmax = 0.1 (horizontal hatches), and a homogeneous distribution with 
Zmax = 1 (\ hatches). For such hard spectra, we need to adjust the normaliza- 
tion energy, Enorm, such that the flux at energies below 10^° eV is consistent 
with observations in this energy range. Thus, we choose the normalization en- 
ergy, Enorm = 4 X 10^^ cV, where the number of events observed by AGASA 
[24] (49) is still statistically signiflcant. 

Again, Figure 7 shows that the number of events aX E > 10^'^ eV is affected 
by the local distribution of the sources. For the adopted normalization, the 
homogeneous distribution gives 3.3±1.6 events above 10^° eV (to be compared 
with the observed 8) and the PSCz distribution provides 3.7 ± 2.0 events in 
the same range. In this last case, about 5% of our realizations give a number of 
events above 10^'^ eV which is equal or larger that the observed one (consistent 
with a 2a signiflcance for a Gaussian error distribution). 

The deflcit of events at energies lower than ~ (3 — 4) x 10^^ eV for hard spectra 
is evident in the case of Zmax = 0.1. Since the high redshift sources contribute 
mainly to the low energy flux, we also considered the case where the maximum 
redshift is Zmax = 1- The increase in Zmax moves the deflcit to energies lower 



11 



26.0 


1 1 1 1 1 1 1 1 1 
- 7 = 3 








- 


25.5 


_ N(>10^°eV)=15.3 + /-3.8 (HOM) 








- 




N(>10^°eV)=17.9 + /-3.5 (PSCz) 












N(>10^^eV)=9075 










T " 


- 


25.0 


: 


<> 


<> T 
<> 


- 




, y .. .L«jfe^^^4», 














- 


E 24.5 


_^- ^^— Y~"Y" P^, 


> 


^^ 


k 








- 


3 


: \ 


% 








~ 


, 24,0 


- 


\ 


, 






- 




- 




• 




- 


23.5 


- 




-. 




. 


23.0 


1 . . . . 






>» 


: 



20.0 
Log[E(eV)] 



Fig. 6. Simulated fluxes for the Auger projected statistics of 9075 events above 10^^ 
eV, and 7 = 3, using a homogeneous source distribution (\ hatches) and the PSCz 
distribution (/ hatches). The solid and dashed lines are the results of the analytical 
calculations for the same two cases. The dash-dotted and dash- dot- dot- dotted lines 
trace the mean simulated fluxes for the homogeneous and the PSCz cases. 
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Fig. 7. Simulated fluxes for the AG AS A statistics of 728 events above 10"*^^ eV, and 
7 = 2.1, using a homogeneous source distribution with Zmax = 0.1 (/ hatches), the 
PSCz distribution with z^ax =0.1 (horizontal hatches), and a homogeneous source 
distribution with Zmax = 1 (A hatches). 

than ~ 2 X 10^^ eV, as shown in Figure 7 (the discrepancy is about a factor 
1.6). 



12 



We also considered the effect of a source luminosity evolution. For instance, 
if the luminosity of the sources increases with redshift z, the flux of UHECRs 
at energies below the GZK cutoff will also increase. This is shown in Figure 
8 for the case m = (horizontal hatches), m = 2 {/ hatches), and m = 4 (\ 
hatches). This type of source evolution by itself does not improve significantly 
the agreement between the theoretical prediction and the AGASA data at low 
energies, leaving unaltered the number of events above 10^'^ eV. In fact the 
height of the bump at ~ (3 — 4) x 10^^ eV also increases with large values of 
m, thus reducing the number of events with energy larger than 10^'' eV when 
the number of events is normalized to the integral flux above 4 x 10^^ eV. The 
number of events above 10^'^ eV is 2.8 ± 1.4 for the case m = 2 and 2.5 ± 1.5 
for m = 4 (these numbers are for a homogeneous distribution of the sources). 




Fig. 8. Simulated fluxes for the AGASA statistics of 728 events above 10"*^^ eV, 
and 7 = 2.1, using a homogeneous source distribution with Zmax = 1 o-nd m = 
(horizontal hatches), m = 2 (/ hatches), and m = A (\ hatches). The lines are the 
results of the analytical calculations, as in [6]. 

In trying to bridge the gap between predicted and observed ultra-high energy 
fluxes, we should consider some additional unknowns that may affect the out- 
put spectrum. For instance, it is not clear how the extragalactic component we 
have considered above gets modified by a Galactic component at the highest 
energies. The spectrum of Galactic cosmic rays might extend up to ~ 10^^ 
eV or even to higher energies. If the Galactic component continues to higher 
energies with the same spectrum as the observed one at lower energies (i.e., 
E~^'^), it will improve the agreement with the AGASA observations when 
combined with the extragalactic flux plotted in Figure 7. We do not have suf- 
ficient information about the Galactic contribution at such high energies to 
rule out this possibility. 
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Another aspect of this problem that should be considered is the effect of 
extragalactic magnetic fields. As shown in [25,26], a large scale magnetic field 
with strength ~ 10~^ - 10~^ G can steepen the spectrum of the observed 
cosmic ray particles significantly. This effect would improve the agreement 
with the AGASA data. But as the diffusive limit is reached, the effect might go 
in the opposite direction: for instance, if the effective diffusion coefficient is too 
small, then the propagation time may be larger than the age of the universe for 
low energy particles generated beyond some distance, considerably reducing 
the cosmic ray horizon [27]. In other words a smaller portion of the universe 
would contribute to the observed flux, thus reducing it. The highest energy 
component of the spectrum should not be affected appreciably by diffusive 
propagation, since most of the particles with E > 10^° eV originate in the local 
universe, where propagation is likely to be non- diffusive. A complete answer 
to this question awaits a full propagation code that includes the effects of a 
realistic model of extragalactic magnetic fields. 



5 Conclusion 



We have studied the effects on the spectrum of ultra-high energy cosmic rays of 
sources associated with a realistic model for the galaxy distribution provided 
by the PSCz and CfA2 catalogs. We considered different injection spectra 
and the possibility of luminosity evolution of UHECR sources with redshift. 
We showed that when the galaxy selection functions are properly accounted 
for, the local overdensity is not as large as found by [8]. Thus, the AGASA 
observations are inconsistent with the predicted flux above f 0^° eV for soft 
injection spectra (7^8). 

We confirm that a local overdensity helps bridge the agreement between theory 
and observations but only (and only slightly) for hard injection spectra (7 ~ 
2). In this case, the observations are within 2a of the predicted number of 
events, in agreement with the findings of [9]. Sources of UHECRs having a 
density field following that of galaxies are therefore consistent with present 
data only for hard injection spectra and, in this case, the GZK cutoff is not 
really a cutoff hut a feature in the overall cosmic ray spectrum. 

As future experiments such as the Pierre Auger Project [28], the High Reso- 
lution Fly's Eye [29], the Telescope Array [30], and the EUSO [31] and OWL 
[32] satellite experiments increase the number of events observed above 10^*^ 
eV, a better determination of the shape of the GZK feature will be obtained. 
The GZK feature will contain information both on the injection spectrum (i.e., 
7) as well as the clustering properties of the sources. Given 7, the clustering 
properties, and the angular distribution of arrival directions, a population of 
sources might be identifiable. If associated with active galaxies or some other 
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specific class of astropliysical objects, tlie sliape of the GZK feature can be 
further used to constrain intergalactic magnetic fields [33] as well as more 
exotic pieces of physics, such as violations of Lorentz invariance [34] . 

Alternatively, the gap between observed and predicted flux may widen as more 
data accumulates. This would indicate that astrophysical proton accelerators 
are unlikely sources of UHECRs. The added difficulty in reaching the extreme 
energies in astrophysical sources further justifies the search for alternative ex- 
planations. Future experiments will play a crucial role in settling this long 
standing mystery, through the determination of the cosmic ray fiux, and ar- 
rival direction distribution, and the great discriminator: the composition at 
extremely high energies. 
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